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Abstract 

We propose a computer-assisted approach to studying the effective 
continuum behavior of spatially discrete evolution equations. The advan- 
tage of the approach is that the "coarse model" (the continuum, effective 
equation) need not be explicitly constructed. The method only uses a 
time-integration code for the discrete problem and judicious choices of 
initial data and integration times; our bifurcation computations are based 
on the so-called Recursive Projection Method (RPM) with arc-length con- 
tinuation (Shroff and Keller, 1993). The technique is used to monitor 
features of the genuinely discrete problem such as the pinning of coherent 
structures and its results are compared to quasi-continuum approaches 
such as the ones based on Pade approximations. 

Mathematical Subject Classification. 65P30, 74Q99, 37L60, 37L20, 
39A11. 

1 Introduction 

In contemporary science and engineering modeling many situations arise in 
which the physical system consists of a lattice of discrete interacting units. The 
role of discreteness in modifying the behavior of solutions of continuum non- 
linear PDEs has recently been increasingly appreciated. The relevant physical 
contexts can be quite diverse, ranging from the calcium burst waves in living 
cells [1] to the propagation of action potentials through the tissue of the cardiac 
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cells [2] and from chains of chemical reactions [3] to applications in supercon- 
ductivity and Josephson junctions [4], nonlinear optics and waveguide arrays 
[5], complex electronic materials [6], the dynamics of neuron chains or lattices 
[7, 8] or the local denaturation of the DNA double strand [9] . 

Whether the phenomenon in question is the propagation of an excitation 
wave along a neuron lattice, the electric field envelope in an optical waveguide 
array, or the behavior of a tissue consisting of an array of individual cells, we 
would often like to model the system through a "coarse level" effective contin- 
uum evolution equation that retains the essential features of the actual (dis- 
crete) problem. Typically computational modeling of such systems involves two 
steps: the derivation of effective continuum equations, followed by their analy- 
sis through traditional numerical tools. In this paper we attempt to circumvent 
the derivation of explicit (closed) continuum effective equations, and analyze 
the effective behavior directly. This is accomplished through short, appropri- 
ately initialized simulations of the detailed discrete process, a procedure that 
we call the "coarse time stepper" . These simulations provide estimates of the 
quantities (residuals, action of Jacobians, time derivatives, Frechet derivatives) 
that would be directly evaluated from the effective equation, had such an equa- 
tion been available. The estimated quantities are processed by a higher level 
numerical procedure (in this case, the Recursive Projection Method, RPM, of 
Shroff and Keller [10]) which computes the effective, macroscopic behavior (in 
this case, traveling waves and their coarse bifurcations). A more general dis- 
cussion of the combination of coarse time stepping with continuum numerical 
techniques beyond RPM can be found in [11]. We have recently demonstrated 
such an approach to the computation of the effective behavior (in some sense, 
homogenization) of spatially heterogeneous problems [12]. This paper consti- 
tutes an extension of this idea to spatially discrete problems. 

The paper is organized as follows: we begin with a brief review of the coarse 
time stepper for spatially discrete problems. We then discuss our illustrative 
problem (a front for a discrete reaction-diffusion system) and its properties. A 
description of our implementation of the coarse time stepper for the bifurcation 
analysis of this particular problem is then presented, followed by numerical re- 
sults. We conclude with a discussion of an alternative approach that involves 
the derivation of an explicit effective evolution equation (based on Padc approx- 
imations), and of the scope and applicability of our method. 

2 A Coarse Time Stepper for Discrete Systems 

Consider a discrete system where each unknown is associated with a point on a 
lattice in space. In the discussion here, we consider a one-dimensional regular 
lattice for simplicity. Higher dimensional and/or possibly irregular, lattices can 
be treated in a similar way. We denote the unknowns {'«(<}• with t € Z. and the 
corresponding points {xg}, such that xg — £Ax, where Ax is the lattice spacing. 
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We assume that the system is governed by the ordinary differential equations 
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(1) 



where n > is an integer representing the range of interaction between lattice 
points. We want to describe this discrete system dynamics through a continu- 
ous function v(t, x) that models the "coarse" behavior of the unknowns on the 
lattice: 



in some appropriate sense. We denote v the coarse continuous solution of (1) 
and we assume that n is not large and that there exists an effective, spatially 
continuous evolution equation for v(x, t) of the form 



for some P and integer M. Such an effective equation for v should "average 
over" the detailed discrete structure of the medium; if there are no macroscopic 
variations of the discrete medium, this equation should therefore be translation- 
ally invariant; for the moment, we will confine ourselves to this case. In terms 
of (1), we can express this as: if F does not depend on £, and if v and v are two 
solutions to the effective equation (2) satisfying v(0,x) = v(0,x + s) for all x, 
then v(t, x) — v(t, x + s) for all time t > 0, all x, and all shifts s. 

It is interesting to consider what the result of integrating such an effective 
equation with a particular, continuum initial condition Vq(x), would physically 
mean. There clearly exists an uncertainty in how such a continuum initial con- 
dition would be imparted to (sampled by) the lattice. One way would be to set 
ui(0) — vo(xi), for all i, but we could equally well set ui(0) = vq(xi + s) for any 
s 6 [0, Ax). There exists, therefore, a one-parameter uncertainty parametrized 
by a continuous shift s. Simulations resulting from different lattice samplings 
of the same continuum initial condition could be quite different. This is best 
illustrated by thinking of a single-peaked function as the continuum initial con- 
dition: the peak may lie precisely at a lattice point, or could fall in-between 
lattice points. It is reasonable to consider as a useful effective continuum equa- 
tion one which takes into account all possible shifts of the initial condition within 
a cell; in analogy with our earlier work [12], we would like to analyze an effec- 
tive equation that would describe the expected result — taken over all possible 
shifts — of sampling the initial condition by the lattice. 

We will use the coarse time stepper approach to simulate an effective equa- 
tion like (2). In this setting, we approximate v(t, x) by the coarse time stepper 
solution u(t, x) at discrete times nT, where T is the time horizon of the coarse 
time stepper. Using the terminology of this framework, we do the following 
steps, starting from a continuous initial condition Vo{x) — u(0,x). 

• Lifting. This initial data (x) is "lifted" to an ensemble of N c different 
initial states of (1) by sampling, 



u e (t) « v(t,x e ), 



v t = P(t,v,d x v,...,d™v), 



(2) 



u J e (0) = v (x£+jAs), 



As = Ax/N ( 



Cl 



j = 0,...,N c -l. (3) 
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Coarse initial condition 



Line up copies and restrict (filter) 




Sample to get shifted copies Integrate each copy independently 




Figure 1: The coarse time stepper: Starting from a coarse initial condition Vo(x), 
lift it by sampling to an ensemble of initial data, {uj(0)}, j = 0, . . . , N c — 1, 
for the system and evolve each set for time T. Line up solutions at time T and 
interpolate to get u(x). Finally filter u(x) to get u(T, x), the result of the coarse 
time stepper at t = T. 
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Setting uj = {u J e }, we write this symbolically as 

where are called the lifting operators. In this case they simply sample 
a continuous function. 

• Evolve. Each ensemble of initial data is evolved till time T according to 
the "true dynamics" (1), 

Uj(T) = T TUj (0), j = 0, . . . , N c - 1. (4) 

where %■ is the solution operator of (1) evolving u(t) to u(t + r). This 
step thus generates an ensemble of solutions Uj(T) at time T. 

• Restrict. Via the restriction operator A4, the ensemble of solutions is 
brought back to a continuous function. 

u(T,x) = M{ Uj (T)}, j = 0,...,N c -l. (5) 

To ensure consistency we require that M{fij) = I. The restriction opera- 
tor M is typically defined as follows. The solutions Uj(T) are thought of 
as sample values of a function u such that u(xi + jAs) = u\. The func- 
tion u is recovered by interpolating the sample values and the restriction 
u(x,T) — M{uj(T)} is finally given as a coarse scale filtering of u(x). 

These steps are illustrated in Figure 1. For n > we define u(nT, x) recursively 
by applying the same construction. Hence, 

u(nT, x) = M{T T Hj}u{{n - 1)T, x). (6) 

The hope is that the coarse time stepper solution u(nT,x), at these discrete 
points in time, can be obtained from a closed evolution equation like (2) whose 
solution, v{t, x) (defined for all t), agrees, at least approximately, with the coarse 
solution obtained from the procedure above, at the discrete points in time, 
v(nT,x) ps u(nT,x). We will refer to the procedure as the coarse time stepper. 

In order to approximate v numerically, we must use a finite representation 
of u{nT,x). We let v n = {wfej^lg 1 , be this representation at time t = nT. The 
elements {v 7 ^ } could be nodal values, cell averages or, more generally, coefficients 
for finite elements or other basis functions. Let II be the operator realizing the 
function from the finite representation, (Jiv n ){x) = u(nT,x). We also require 
that the restriction operator projects on the subspace spanned by the finite 
representation, and we can redefine it to also convert the projected function to 
this representation. Symbolically, we then write the coarse time stepping 

v n+i = M {r T ^}Uv n =: G(v n ). (7) 

Note that we may not be able to write down the explicit expression for G or 
the equation (2) for v(t, x), but our definition of u(t, x) allows us to realize its 
time-T map numerically in a straightforward fashion. 
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Applied directly to the simulation, the coarse time-stepper does nothing to 
reduce the cost of detailed computation with the discrete dynamics. It is only in 
conjunction with other techniques (like projective integration [13], or matrix-free 
fixed point techniques) that the coarse time stepper may provide computational 
or analytical benefits. Here we will make use of the coarse time stepper in 
conjunction with the Recursive Projection Method (RPM), to perform stability 
and bifurcation analysis of certain types of solutions of the (unavailable) coarse 
evolution equation. For a schematic illustration of the coarse time stepper with 
RPM, see Figure 2. 
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Figure 2: An overview of the coarse time stepper with RPM. 

RPM helps locate fixed points, allows us to trace fixed point branches and 
locate their local bifurcations; when the bifurcations in (7) that we are interested 
in do not involve fixed points, G has to be reformulated. How this is done 
depends on the application; for the type of solutions considered here (traveling 
fronts), the appropriate modification is discussed in Section 3.2. 



G 



3 A Discrete Traveling Front Example 



The effects of discreteness on the propagation of traveling wave solutions have 
been documented and analyzed in many different settings over the last two 
decades. From the pinning of travelling waves in discrete arrays of coupled 
torsion pcndula and Hamiltonian models [14, 15], to the trapping of coherent 
structures in dissipative lattices of coupled cells [16, 17] (see also references 
therein) , the role of spatial discreteness has triggered a large interest in a diverse 
host of settings. 

Effective equations capable of describing the nature of the solutions of dis- 
crete problems should successfully capture the effects of discreteness on the 
traveling wave shape and speed. More importantly, they should be capable of 
accurately predicting qualitative transitions (bifurcations) that are inherently 
due to the discreteness. The most prominent of those is probably the pinning 
of traveling waves and fronts often observed when the lattice spacing becomes 
sufficiently large. To illustrate the performance of our proposed coarse equation 
in capturing such a front pinning, we have chosen what is arguably a prototypi- 
cal spatially discrete problem capable of exhibiting it: a one-dimensional lattice 
with scalar bistable on-site kinetics and nearest neighbor diffusive coupling be- 
tween lattice sites. Our test problem is, therefore, a discrete reaction-diffusion 
system described by 

^ = ^2(«/-i-2«/ + «/+i) + /(u/), ieZ (8) 

with 

f(u) = 2u(u-l)(ri-u), 77 = 0.45. (9) 

This can serve as a model of e.g. individual cells in the cardiac tissue 
which are resistively coupled through gap junctions (see e.g., [16] and refer- 
ences therein). In this case the solution U£, would correspond to the electrical 
potential of the cells. For small Ax the system possesses solutions that can be 
characterized as discrete traveling fronts: see Figure 4. These solutions have 
a near constant shape and travel in a "lurching" manner. When Ax becomes 
sufficiently large, front propagation fails (front pinning). In our example, this 
happens at Ax — Ax* s=s 2.3, see Figure 3. The front speed for an infinite lattice 
approaches the asymptotic "PDE speed" value 0.1 as the lattice size tends to 
zero. 

We will examine how faithful the coarse time stepper is to the properties 
of the solutions of the full discrete model (8). Our numerical simulations are 
restricted to a finite domain, using N — 64 grid points. At the boundaries, we 
prescribe Neumann-type conditions 

UN - U N -i = 0, 

uq — U-i = 0. 

This should model the full problem accurately as long as the (relatively narrow) 
front is positioned sufficiently far from the boundary. 
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Speed versus grid partioning 




Figure 3: The speed of the front as a function of Ax. As the lattice spacing is 
increased, the speed v approaches zero; the front stops at Ax* w 2.3. 
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(a) u(t, x) in the a;t-plane 



(b) u(t,x), t = 0, 2.5, 5, . . . , 40 



Figure 4: The plot illustrates how the front advances when Ax = 1.75. The 
left figure shows the front in the xt-plane; the grayscale is proportional to the 
solution u(t, x). The right figure shows the solution as a function of x at different 
time levels. The time interval is t E [0,40]. Looking at the spacing between the 
solution instances, we can see how the front speed varies in a lurching manner. 
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3.1 Construction of the coarse time stepper 

In this section we detail the procedures associated with the coarse time stepper 
applied to the test problem (8, 9) on the finite interval I = [0, L], where L = 
NAx and the cell locations are Xj = jAx, with j = 0, . . . , N — 1. 

Our choice of finite representation of the coarse solution are M nodal values 
v n = k = 0, ...,M — 1, evaluated at t = nT and y^ = fcAy, with 

MAy = NAx. 

For many solution shapes Fourier interpolation would be a natural interpo- 
lation operator realizing the coarse solution u(nT, x) from v n . We denote direct 
Fourier interpolation by 11^. We could then define the corresponding lifting 
operators \J- via the shifting operator <S/ : M M — > R N , 

^u:=S } ]As u, (S{u) e := (U f u)(x e + s), s > 0, 

where IT^ uses {yt} as interpolation nodes. In our case, however, the solution 
is not periodic on / and we get large errors if we use 5/ directly. Instead we 
apply Fourier interpolation to the differences of the v n sequence. We thus use 
the modified shifting operator S s : R M — > R N given by 

e 

S s u := CS{D, (Cu) e := 1 + ]T Uj , {Du) e := 

3=0 

(10) 

We then define the lifting operator \i : M M — > M. NxNc (acting directly on v n ) as 
l» B = {ft« B }, N v n :=S jAs v n , As = ^, 

where j = 0, . . . ,N C — 1. 

The restriction operator A4 : M. NxNc — > R M is also defined using the shifting 
operators, but now with negative shifts, 

S f _ s :R N ^R M , (S f _ s u) k :=(Uf U )(y k -s), s>0, 
where ±1^ uses {xi} as interpolation nodes. We then set S- s = CSL S D and let 

j iVc-l 

Note that these choices of (i and M are consistent when TV > M. Then, by 
the sampling theorem SL S S? = I on R M . Moreover, it is easy to see that 
CD = DC = I. Therefore, we also have 

S- S S S = CSL s DCS f s D = CSL a S J a D = CD = I, 

on R M and consequently, 

N c -1 N c -1 N c -1 

M ^ n = AT E <5-,a sW «" = W E 5-^,5^^" = — ]T »" = 

C 3=0 C 3=0 C 3=0 



uo-1, £ = 0, 

£ > 0. 
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We should also remark here that, in the special case when N = M, we have 



(P N U f u)(xe), u = {u r }, U£ +jN = it 



j 



where Pn is a projection on the N lowest Fourier modes. Hence, if we used direct 
Fourier interpolation and M — N, then our definition of M is equivalent to 
lowpass filtering of it, the lined up copies described in Figure 1, top right. When 
we replace <S/ by S s we do not retain exactly this property, and a definition of M. 
based on simple lowpass filtering is no longer consistent. However, our procedure 
still corresponds to a type of lowpass filtering, although a more complicated one. 

For the time integration of (8) we use the Crank-Nicolson method, treating 
the nonlinear term explicitly Thus, with u; = {w®} 6 R N , 

T T w° := w Nt = {wf T }, N T At = T, 

where {w"} are given iteratively by 

for I = 0, . . . , N — 1, together with the free boundary conditions 

w^-wZ = 0, 
w% ~ w , x_ 1 = 0. 

In our computations we use the time step At = 0.01. 



3.2 Steady state formulation 

The coarse solution u(nT,x) as we have defined it is a (practically) constant 
shape moving front. In order to convert this moving state into a stationary state, 
we can factor out the movement through a procedure based on template fitting 
([18, 12], see also [19]) which pins the traveling front at a fixed ^-coordinate. 
This is performed by a "pinning-shift" operator, which we denote V . Our coarse 
time stepping is then modified from (7) to 

v n+i = VM{r T ^}Uv n =: G(v n ). (11) 

This formulation has a steady state at the constant shape moving front. 

Let us start from the basic, Fourier based, pinning-shift operator V* : M M — ► 
M M . After introducing a template function S(x), we define 

V s w := S{ w, c — argmax / (U f w)(x + c r )S(x)dx. (12) 

c'eR Jo 
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Hence, P?w is the shifted version of w that best fits the template S(x), in the 
sense that it maximizes the L2-hiner product between its Fourier interpolant 
and S. Upon convergence, the effective front speed v can be deduced from the 
converged value of c and the time reporting horizon T simply by taking v = c/T . 
With the template S(x) = 1 — cos(2irx / L) we can compute the inner product 
in (12) explictly, 

^■J (n f w)(x + c r )S(x)dx = w Q -3?(w ie lc ') , (13) 

where Wk are the Fourier coefficients of w. Hence, since wo is real c in (12) 
should be chosen such that w\e rc is real and negative. This is easily implemented 
numerically together with the Fourier shift S[. 

For the same reasons as in the implementation of the coarse time stepper, 
we would like to avoid direct Fourier interpolation of the solution, since it is 
not periodic. Therefore, we modify V* to operate on differences instead. In the 
same spirit as in Section 3.1, we let 

V := CP f D, 

with C and D defined in (10). We still use the effective propagation speed given 
by pf. 

An important property of the Fourier based pinning shift operator is that it 
satisfies {P^) 2 — V* , which follows from the sampling theorem [12]. For other 
types of interpolation, such as piecewise polynomial interpolation, the pinning 
shift operator will not have this property and a steady moving coarse shape may 
not translate into a fixed point for (11). Our modification still has this property 
though, since 

P 2 = CP f DCP f D = C(P f ) 2 D = CP f D = P, 
where we used the fact that DC = I. 



3.3 The RPM with pseudo-arclength continuation 

RPM is an iterative procedure which can accelerate the location of fixed points 
of processes; under certain conditions it can help locate steady states of dynamic 
processes (in particular, discretized parabolic PDEs). It can be an acceleration 
technique for the solution of nonlinear equations, and a stabilizer of unstable 
numerical procedures (as it was first presented, [10]). Consider the fixed point 
problem 

F(u;X) = u, (14) 

and let J be the Jacobian of F. 

• Like the Newton method, RPM can converge rapidly to the fixed point so- 
lution u* provided the initial guess is good enough; the convergence occurs 
even if J(u*) has a few eigenvalues larger than one. The computational 
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cost and convergence rate depend on the eigenvalues of J. Optimally there 
should be a clear gap in the spectrum between small and large (near the 
unit circle) eigenvalues and a limited number of large (in norm) eigenvalues 
for RPM to perform well. 

• J never needs to be evaluated directly, only F. We can therefore apply 
RPM to any "black box" code that defines a function F; it is a "matrix- 
free" method. 

• As a by-product, RPM also computes approximations of the largest eigen- 
values of J. This gives approximate stability information about the fixed 
point. 

When RPM is used for the computer-assisted bifurcation analysis of steady 
states of (usually dissipative evolution) PDEs, the function F represents a time- 
stepper, a subroutine that takes initial data and reports the solution of the PDE 
after some fixed time (the reporting horizon T). A fixed point then satisfies (14). 
The conventional way of finding the steady state using a time-stepper would be 
to call it many times in succession — in effect, to integrate the PDE for a long 
time, corresponding to solving (14) by simple fixed point (Picard) iteration. 

RPM can improve this approach in two important respects. First, the con- 
vergence can be significantly accelerated. The nature of many transport PDEs 
usually encountered in engineering modeling (the action of viscosity, heat con- 
duction, diffusion, and the resulting spectra) dictates that there exists a separa- 
tion of time-scales, which translates into an eigenvalue gap in the spectrum of J 
at the steady state. Second, RPM converges even if the steady state is slightly 
unstable, i.e. when J has a few eigenvalues outside the unit circle. It may thus 
be possible to compute (mildly) unsteady branches of the bifurcation diagram 
using forward integration (but in a non-conventional way, dictated by the RPM 
protocol). RPM still retains the simplicity of the fixed point iteration, in the 
sense that no more information is needed than just the time-integration code. 
This code, which may be a legacy code, and can incorporate the best physics 
and modeling available for the process, is used by RPM as a black box. 

RPM can be seen as a modified version of fixed point iteration. It adap- 
tively identifies the subspace corresponding to large (in norm) eigenvalues of 
J, hence the directions of slow or unstable time-evolution in phase space. In 
these directions the fixed point iteration is replaced by (approximate) Newton 
iteration. More precisely, suppose F : R N x M — > R N in (14). Let P be the 
maximal invariant subspace of J corresponding to the m largest eigenvalues and 
let Q be its orthogonal complement in R N . The solution u is decomposed as 
u = p + q = Pu + Qu, where P and Q, are the projection operators in R N on 
P and Q. These are constructed from an orthogonal basis V p 

P = v p v*, 

Q = I-V P V?. 

In a pseudo-arclength continuation context the solution u = u(s) and A = A(s), 
where s parameterizes the bifurcation curve. In addition to (14) we then use an 
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algebraic equation to be able to handle turning points, 



S(u, A, As) 



\u(s) u(s Asjf + |A(s)-A(s-As)| 2 _ ^ 



As 



As 



0, (15) 



where u(s — As) and A(s — As) refers to the converged solution at the previous 
point on the continuation curve. 

The solution is advanced using a predictor-corrector method. Via extrap- 
olation from previous points Ui = u(si), Xi — A(sj) and As^ = Sj+i — s,, the 
predictor-solution is obtained. Comparing a first order extrapolation, 



A* = Xi + 



with a second order extrapolation, 



A; — \{- 



Asi_i 
Ui - Ui-i 



-A Si 



As, 



As i5 



l A,(l- 7 )-2A t _ 1 + (l+ 7 )A t _ 2 A 2 
2 Asi_iAsj_ 2 

1 Uj(l - 7) - 2ll i _ 1 + (1 + 7)Ui_2 



and requiring that 



u + „ 

2 Asi_iAsi_ 2 

Asj-i - Asj-2 

As^i + Asi_2 ' 

max(||u** -u*||,|A** - A*|) < e 



As? 



(16) 



the stepsize is determined. Here e is a user specified tolerance. As the corrector 
method, we use RPM with pseudo-arclength continuation, see [10, 30]. Starting 
from u° = u** and A = A**, the iterative scheme is given by 



slv p 



I) V p T F x - 




Ap ' 




Sx 




AA 





u n+1 

X n+1 



QF(u n ,\ n ), 

V p T F(p n 

S(p n + q 

p n + V p Ap n + q n+1 
A" + AA", 



A") 

n+1 \n\ 



,A") 



where the left hand side consists of partial derivatives of S in (15) and of F in 
(14) with respect to u and A. The iterates u n — p n + q n will converge to the 
solution of (14) under the assumptions discussed above. If the number of large 
norm eigenvalues, to, is limited, the dimension of F and the projected Jacobian 
in the Newton iteration, Vj JV P — I, remains small. Only this small matrix 
needs to be inverted. For a more complete description of RPM we refer to [10]. 
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4 Numerical Results 



In this section we present some numerical results using the coarse time stepper 
and the procedure described above to simulate an effective equation for the 
discrete problem in (8). We will start by discussing the "exact" bifurcation 
diagram of the discrete system, which we attempt to approximate. We will then 
show results obtained through the coarse time stepper, and discuss the effect of 
time stepper "construction parameters" like the reporting time horizon, T (the 
time to which (1) is integrated within the coarse time stepper), and the number 
of different initial shifted copies, N c . 

Figure 5 shows the bifurcation diagram of the discrete problem as a function 
of the parameter Ax, the lattice spacing, in the regime close to the onset of 
pinning. For lattice spacings smaller than Ax* « 2.3 the system has, as we 
discussed, an attracting, front-like solution that travels; its motion is modulated 
as it "passes over" the lattice points. For an infinite lattice, this modulated 
traveling solution possesses a discrete translational invariance: ue + i(t + r) = 
Ui(t). The shape of the modulating front is shifted by one (resp. 2,3, ... ,n) 
lattice spacing after time r (resp. 2r, 3r, . . . , nr); this helps us define its effective 
speed v(Ax) = ^ (see Figure 3). As Ax approaches zero, for an infinite lattice, 
the discrete front approaches the continuum front of the PDE, and its speed (the 
period of the modulation divided by Ax approaches the PDE front speed, 0.1 
(see see Figure 3). 

If we identify shapes shifted by one lattice constant, the attractor appears as 
a limit cycle with period r. As the lattice spacing approaches the critical value 
Ax* the speed of propagation approaches zero (the period of the "limit cycle" 
approaches infinity); asymptotically, v(Ax) ~ \Ax — Ax*|°- 5 . As discussed in 
[20, 21] what occurs is a Saddle-Node Infinite Period (SNIPER) bifurcation: 
a saddle-node bifurcation where both new fixed points appear "on" the limit 
cycle. For larger values of Ax the "saddle" and the "node" move away from each 
other, and what used to be the limit cycle is now comprised from the saddle, the 
node, and both sides of the one-dimensional unstable manifold of the saddle, 
which asymptotically approach the node. 

The saddle and the node are, of course, stationary fronts. A pair of them 
exists for every "unit cell": all "node fronts" are shifts of each other by one 
lattice spacing, and all "saddle fronts" are also shifts of each other by one lattice 
spacing. Since the medium has a discrete translational invariance, this makes 
sense — if an initial condition gives rise to a front eventually pinned at some 
location in the discrete medium, the shift of this initial condition by one lattice 
spacing will eventually get trapped one lattice spacing further. This saddle- 
node bifurcation can be seen in Figure 5a; linearizing around the saddle front 
will give a positive eigenvalue X s , while the corresponding eigenvalue A„ for the 
node front would be negative. Since we look at the problem in discrete time, 
what is plotted is the multiplier ^t„ ;S = cxp(A„, s T), where T is the reporting 
horizon. The saddle front has a multiplier larger than 1, while the corresponding 
multiplier for the stable node is less than 1; both multipliers asymptote to 1 at 
the SNIPER (Ax*). 
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Figure 5b shows the bifurcation diagram in terms of the front traveling speed. 
Since both the saddle and the node fronts are pinned (have zero speed) they 
both fall on the zero axis; we plotted their eigenvalues in Figure 5a to distinguish 
between them. The true traveling speed (broken line) is compared with the 
effective traveling speed predicted by a coarse time-stepper using N c = 5 copies 
within each unit cell, and a reporting horizon of T = 32. The coarse time stepper 
speed is a byproduct of fixed point computation and continuation with it; short 
bursts of detailed simulation arc used in the RPM framework to construct a 
contraction mapping that converges to a fixed point of the time stepper. The 
final shift upon convergence (from the pinning-shift computation), divided by 
the time stepper reporting horizon gives us an estimate of the "effective speed" . 
Inspection of Figure 5b indicates that the coarse time stepper never predicts 
a speed that is exactly zero; yet it gives a good approximation of the effective 
speed, all the way from small Ax to the near neighborhood of the pinning 
transition, when the effective speed becomes small. 

We will return to discussing this issue of "small residual motion" for the 
coarse time stepper shortly. To give an indication of when the procedure stops 
being quantitative, we have included the Ax/T curve in Figure 5b: disagreement 
starts well in the regime where the effective movement is less than one unit cell 
per observation period. In the next section we will compare the "goodness of 
approximation" of our coarse time stepper to the effective speed predicted by 
the Pade approach to extracting effective continuum equations. It is interesting 
that the coarse time stepper sometimes predicts a small hysteresis loop at low 
speeds, relatively close to "true pinning" ; notice in Figure 5a the unstable (larger 
than one) multipliers for the brief saddle part of this loop. We will discuss a 
tentative rationalization of this below. 

Figure 6 illustrates the effects of "time stepper construction" parameters on 
the effective behavior predicted by the time stepper: the reporting time-horizon, 
for two different sets of shifted copies (N c = 3 and N c = 10) as well as the effect 
of the number of copies for a fixed time horizon (T = 16). Augmenting the 
time stepper reporting horizon is shown in Figure 6a-b; clearly, in both cases, 
extending the time stepper reporting horizon extends the region over which its 
effective speed agrees with the true problem closer to Ax*. Larger numbers 
of copies (N c — 5, 10, 20) also perform slightly better than smaller numbers 
(N c = 3). In all cases the qualitative behavior is the same: (a) successful 
approximation of the effective speed until reasonably close to true pinning; (b) 
all differences occur when the average front motion is significantly less than 
one unit cell per reporting horizon; (c) there is always a slight residual motion, 
which — possibly after a small hysteresis loop close to true pinning — eventually 
becomes negligible. 

We now turn to the discussion of the slight residual motion of the coarse time 
stepper at large Ax beyond Ax*. For an infinite domain, the saddle and node 
pinned fronts appearing there are invariant to translations by one lattice spacing; 
for a large enough computational domain we still see two pinned front solutions 
per cell. When we "sprinkle" initial conditions along the cell, depending on their 
location with respect to the saddle front, the trajectories may either be attracted 
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to the stable node "to the right" or to the one "to the left" of the saddle. It is 
instructive to represent these solutions as in Figure 7a, in a way that identifies 
the "right" node front with the "left" one; here translation along the lattice 
corresponds roughly to rotation along the circle. The node is denoted by a black 
circle, and the saddle by a white one. The small squares represent the initial 
positions of our initial condition "copies" . The fate of our distribution of initial 
conditions is governed by their initial "angle" on the circle — as our time horizon 
grows all initial conditions will asymptote to a stable front, cither the left one 
(moving counterclockwise on the circle) or the right one (clockwise movement). 
We now see clearly the physical reason behind the net residual motion for any 
finite time horizon for the coarse time stepper. An initial condition that is put 
down "at random" in a unit cell deep in the pinned regime, even if it never exits 
this unit cell, will gradually traverse the part of the circle separating it from the 
closest node front. 

When the critical parameter value is approached from the pinned side, the 
saddle and the node fronts approach each other on the circle, on their way to 
coalescing at the SNIPER bifurcation point [20, 21). Figure 7b shows how this 
process becomes manifest in the coarse time stepper computations, using the 
problem in Figure 5 as our example. Deep in the pinning regime (high Ax, 
marked a) the relative "phase" of the saddle and the node pinned fronts on the 
circle remains roughly constant. The distance each member of our ensemble of 
initial conditions has traversed during one time horizon can be deduced from 
Figure 7b: the copy with the largest negative movement is the one closest to the 
saddle but on its left (copy number two). One can similarly rationalize the la- 
belling of the remaining curves in Figure 7b. When Ax is reduced approaching 
the onset of pinning, at some point the saddle front starts moving apprecia- 
bly towards the node front. As part of this movement, it "sweeps" the circle 
counterclockwise; at Ax « 2.8 it has its first encounter with one of our initial 
conditions — the closest one on the left. When the saddle "moves past" it into 
the regime marked (3, this copy, which was responsible for the largest negative 
displacement now approaches asymptotically the node front on the right, per- 
forming the largest positive displacement (and so on for the remaining copies). 
Eventually, in the propagating regime, marked 7, and for long enough reporting 
horizons, the initial "phase" difference (a fraction of a cell) becomes negligible 
compared to the net displacement of each point (several cells). 

The real movement in phase space is shown un Figure 7c for two different Ax. 
In these subfigures, the x-axis represents sin(27ra; c ) where x c corresponds to the 
location of the front, more specifically x c — ^2 e (Du)i£. The y-axis represents 
max; \ (Du)g\. The initial positions of the copies are indicated by small squares 
and their locations at t = T, the time horizon, are marked by filled circles. The 
labels refer to the same copies as in Figure 7b. 

As the reporting time horizon of the time stepper goes to infinity, it is clear 
that one can compute the average residual movement from the asymptotic posi- 
tion of the saddle front, i.e. from the relative extent of the circle "to the right" 
and "to the left" of the saddle front. The most reasonable point to "declare" as 
an estimate of the true pinning from coarse time-stepper computations would 
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come from a polynomial extrapolation of the "successful" regime (close to the tip 
of the "apparent parabola" in Figure 5); alternatively, a value of Ax where the 
speed is small enough (well below one unit cell per time horizon) and its varia- 
tion with number of copies and time horizon is below a user-prescribed tolerance, 
would also serve this purpose. While there is no well defined pinning bifurcation 
for the coarse time stepper (since pinning is an inherently non-translationally 
invariant bifurcation), the procedure can provide a good approximation of the 
effective shape and speed of the traveling fronts, as well as "common sense" 
ways of numerically estimating the true pinning. 

5 An Alternative Continuum Approach: Pade 
Approximations 

In this section, we propose an alternative scheme for capturing effects of discrete- 
ness, by means of a (now explicit) continuum equation. This PDE is obtained 
by means of Pade approximations [22, 23] which can be used to approximate 
discreteness in a quasi-continuum way, through the use of pseudo-differential op- 
erators. In particular, starting from the Taylor expansion for analytic functions, 
see e.g., [24], 

u(x + m) = exp(md x )u(x) , 
one can then express spatial discreteness as 

ii£ + i + u£-i — 2u£ = (exp(Axd x ) + exp(— Axd x ) — 2)u(x) 
= 4sinh 2 (^-) U (x,i). 

Expanding exp(±Axd x ) [23], one then obtains 

1 At 2 1 

exp(±Axd x ) - 1 = -Ax 2 (l + —d 2 x + ...)dl± Ax{l + -Ax 2 d 2 x + ...)d x 

Finally, regrouping the terms in the manner of Pade [22, 23] yields 

, , a „ x , 1 A x 2 dl Axd x 
exp(±Axd x ) - 1 « ± (17) 

We now use the pseudo-differential operator approximation in (17) to convert 
the discrete model in (8) into the PDE approximation of the form: 

12 x 

Such approaches were introduced and used extensively by Rosenau and collab- 
orators [25, 26, 27] to regularize nonlinear wave equations, particularly of the 
Klein-Gordon type. 
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Figure 7: Movement of the individual copies, for N c = 5, T = 32. 
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Eq. (18) clearly emulates the discrete setting in some key aspects of the rel- 
evant spectral operator properties (i.e., of the discrete Laplacian in comparison 
with the pseudo-differential operator of (18)). For example, considering plane 
wave solutions of the form exp(Ai — ikx), wc obtain in the discrete case the 
linearized dispersion relation (around a uniform state u — Uhom) 

2 

A = ^ (cOs(kAx) - 1)) + /'(Uhom)- 

In the case of (18), the corresponding equation becomes 

k 2 



l + ^k 2 



+ f'(u ho 



Apart from sharing the continuum limit, the two dispersion relations share 
another qualitative feature which is particularly important [25, 26, 27]; namely, 
the presence of a lower bound in the continuous spectrum. Notice, however, 
that the two lower bounds are different (/'(uhom) — 4/ Ax 2 in the discrete case 
versus /'(whom) — 12/ Ax 2 in the Pade approximation). 

It would then be of interest to alleviate this spectral discrepancy, as well 
as to match the discrete operator (if possible) to a higher order in the Taylor 
expansion 

u n+ i + un-i -2u n _f, 2Ax 2 ^ 2 {2j) Ax^ + 0(Ax e ) 

Ax 2 ~ 4-*> (2j)\ ~ U xx+ 12 U xxxx -f U Sx -f U{L^X j. 

This can be achieved by a natural generalization in the form of a continued 
fraction such as e.g., 

d 2 

(19) 



In order to use (19) in practice (i.e., for computational purposes), we convert 
the three fractions into one of the form 

d 2 x (l + aAx 2 d 2 ) 



1 + (a + (3)Ax 2 d 2 + -/Ax^f 



where a simple (algebraic) reduction of A,B,C to a, (3, 7 has been used. We 
then use Taylor expansion of the denominator to convert the expression of (20) 
into one resembling (19). By matching up to 0(h 6 ) the exact Taylor expansion, 
we obtain three algebraic equations for a, (3 and 7. In this way, we obtain a set 
of solutions for a, (3 and 7. We use here the set a — —0.007 912, (3 = —1/12, 
7 = 0.002 056. An additional benefit (to the matching of the Taylor expansion 
up to correction terms of 0(h 8 )) that should be highlighted here is the value 
a/(-iAx 2 ) = 3.848/Ax 2 of the lower bound expression for A, which is much 
closer to the theoretical lower bound of 4/ Ax 2 than the prediction 12 /Ax 2 of 
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the leading order approximation presented previously. The resulting evolution 
equation will then read: 

dZ(l + aAa?dl) 
Ut ~ 1 + (a + /?) Ax 2 ^ + 7 Ax 4 94 + 1 W W 

Both (18) and (21) can be numerically implemented in a straightforward 
manner, by means of the spectral techniques described in [28]. We have per- 
formed numerical simulations of the front propagation, using 1024 modes in the 
spectral decomposition of (18) and (21). We will refer to these equations as the 
(Pade) models A and B respectively. A fourth order Runge-Kutta algorithm 
has been used for the time integration. For each value of Ax, we identify the 
position x c of the front as the point where the ordinate of the front acquires 
the value u = 1/2. The linear interpolation scheme suggested in [29] has been 
implemented and has proved to be an efficient front tracking algorithm in all 
the examined cases. 




1 1.5 2 2.5 3 3.5 4 

Ax 

7 

Figure 8: Effective front speed as a function of Ax, for the Pade model A (top 
panel) and model B (bottom panel). 

Our results of this quasi-continuum approach to the discrete problem can 
be summarized in Figure 8 and Figure 9. Figure 8 shows the speed of the 
fronts in Pade models A and B respectively. We can observe that the critical 
value of Ax beyond which trapping of the front occurs is significantly displaced 
from the actual one of Ax* « 2.3, for r\ = 0.45. In particular, for model A, 
Ax* « 6.4, while for model B, the corresponding critical value is Ax* « 3.8. 
We can deduce that the latter model is closer to the actual physical reality, even 
though the relevant prediction is still considerably higher than its actual value 
for the discrete model. 
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X/AX 

Figure 9: The figure shows the time evolution of the front for model A and 
for Air = 6.4. The top panel shows the time evolution of the front center 
which eventually leads to trapping. The bottom panel shows the final front 
configuration of the numerical simulation at t = 150. 

In part at least, these results (and the discrepancy from the actual discrete 
case) can be justified by observing Figure 9. The bottom panel of the figure 
suggests that the only way in which the front can stop in these quasi-continuum 
Pade approximations is by becoming practically a vertical shock-like structure. 
In this case, the "mass" of the front which is given by u^dx (see e.g., [29] 
and references therein) becomes practically infinite. This means that the inertia 
of the front becomes too big for the front to move and hence "pinning" occurs. 
However, notice that this process of pinning is significantly different than the 
details of the discrete structure of the problem (such as e.g., the saddle- node 
bifurcation and the transition to pinned solutions) . The translationally invariant 
quasi-continuum Pade approximations of models A and B do not "see" such 
features. Instead, they incorporate the well-known feature of front steepening 
for stronger discreteness [15] and the criticality of the latter feature eventually 
leads to pinning. 

An additional pointer to the fact that such (pseudo-differential operator) 
models are "eligible" to pinning is that they are devoid of some of the impor- 
tant symmetries that are inherently related to traveling such as the Galilean 
invariance in the case of continuum bistable equation or the Lorentz invariance 
of its Hamiltonian (nonlinear Klein-Gordon) analog. 
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6 Summary and Discussion 



We presented a computer-assisted approach for the solution of effective, trans- 
lationally invariant equations for spatially discrete problems without deriving 
these equations in closed form. Assuming that such an equation exists, its time- 
one map is approximated through the coarse time stepper, constructed through 
an ensemble of appropriately initialized simulations of the detailed discrete prob- 
lem. Combining the coarse time stepper with matrix-free based numerical anal- 
ysis techniques, e.g. contraction mappings such as RPM, can then help analyze 
the unavailable effective equation. We are currently exploring the use of our 
coarse time stepper with coarse projective integration [13, 11, 31]. Matrix-free 
eigenanalysis techniques should also be explored, especially since they can help 
test the "fast slaving" hypothesis underlying the existence of a closed effective 
equation (see, for example, the discussion in [32, 33]). 

We also presented initial computational results exploring the effect of certain 
"construction parameters" of the approach: the number of shifted copies in the 
ensemble of initial conditions, as well as the time-horizon used. We included 
a comparison between our approach and a particular way of obtaining explicit 
approximate translationally invariant evolution equations for such a problem 
(the Pade approximation). More work is necessary along these lines, exploring 
the relation of our approach with traditional homogenization methods at small 
lattice spacings. A discrete problem whose detailed solution can be obtained 
explicitly (perhaps a piecewise linear kinetics problem) or at least approximated 
very well analytically over short times, would be the ideal context in which to 
study these issues. 

Several extensions of the approach can be envisioned, and might be inter- 
esting to explore. A time stepper based approach can be applied without mod- 
ification to hybrid discrete-continuum media, e.g. continuum transport with a 
lattice of sources or sinks, such as cells secreting ligands into and binding them 
back from a liquid solution, [34]. It is clear that it can be tried in more than 
one dimensions, and for regular lattices of different geometry. For irregular lat- 
tices the averaging "over all shifts" we performed here for periodic media can 
be substituted with a Monte Carlo sampling over the distribution of possible 
lattices that takes into account what we know about the statistical geometry of 
the lattices. In this paper we assumed that an equation existed and closed for 
the expected shape of the solution. Conceivably one can attempt to develop time 
steppers not only for the expectation (the first moment of a distribution of possi- 
ble results), but, say, for the expectation and the standard deviation of possible 
results; the lifting operator would then have to be appropriately modified. Fi- 
nally, our time stepper here was built on short simulations of the entire detailed 
discete system in space. Hybrid simulations, where a known, explicit effective 
equation is accurate over part of the physical domain can be done; an "overall 
hybrid coarse" time stepper (explicit equation over part of the domain, and the 
coarse time stepper in this paper over the rest of the domain) will then be used. 
In a multiscale context, we have proposed "gaptooth" and "patch dynamics" 
simulations [35, 36], where the present coarse time stepper integrations are per- 
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formed not over the entire domain, but over a mesh of small computational 
"boxes". Both hybrid and "gaptooth" simulations, if possible, require careful 
boundary conditions for the "handshaking" between the continuum equation 
and the discrete simulations, or the discrete simulations in distant boxes, ef- 
fectively implementing smoothness of the solution of the unavailable effective 
equation (e.g [36, 37, 38, 39, 40]). 

We close with a discussion of the "onset of pinning" , the transition around 
which our test example of the coarse time stepper was focused. Continuum 
effective equations such as the ones discussed here through the numerical time- 
stepping procedure do not, strictly speaking, possess a bifurcation at the critical 
point of the genuinely discrete problem. In this effective process, the bifurca- 
tion is smeared out and rendered a "continuum transition" (see, for example, 
materials science models of the onset of movement of a front, [41, 42]). On the 
other hand, one might argue that this is an acceptable, and possibly optimal 
way for a continuum equation to represent the discrete bifurcation to pinning. 
We can see that other procedures, such as the discreteness-emulating Pade type 
ones, lose a lot of the quantitative structure of the relevant transition. On the 
other hand, if a continuum differential (as opposed to pseudo-differential) equa- 
tion was constructed to "model" this transition, the latter would possess other 
artificial features such as a topologically mandated, unstable branch of traveling 
wave solutions [43]. It is conceivable that the short hysteresis loop sometimes 
predicted by the coarse time stepper close to pinning conditions is a "vestige" of 
this unstable branch that translationally invariant equations would necessarily 
predict. In conclusion, it can be appreciated that genuinely discrete problems 
and continuum ones have inherent differences 1 that cannot be fully captured by 
emulating (or "summarizing") the one context through the other. Nevertheless, 
the approach proposed here, combined with a "common sense" interpretation 
of its results with respect to the genuinely discrete problem, performs in a sat- 
isfactory way for the modeler, even for the "most different" features between 
discrete and continuum models. 
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